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Abstract The recently developed essentially fourth-order or higher low dissipative 
shock-capturing scheme of Yee, Sandham and Djomehri (1999) aimed at minimizing nu- 
merical dissipations for high speed compressible viscous flows containing shocks, shears 
and turbulence. To detect non smooth behavior and control the amount of numerical 
dissipation to be added, Yee et al. employed an artificial compression method (ACM) 
of Harten (1978) but utilize it in an entirely different context than Harten originally 
intended. The ACM sensor consists of two tuning parameters and is highly physical 
problem dependent. To minimize the tuning of parameters and physical problem depen- 
dence, new sensors with improved detection properties are proposed. The new sensors 
are derived from utilizing appropriate non-orthogonal wavelet basis functions and they 
can be used to completely switch off the extra numerical dissipation outside shock lay- 
ers. The non-dissipative spatial base scheme of arbitrarily high order of accuracy can 
be maintained without compromising its stability at all parts of the domain where the 
solution is smooth. Two types of redundant non-orthogonal wavelet basis functions 
are considered. One is the B-spline wavelet (Mallat & Zhong 1992 ) used by Gerritsen 
and Olsson (1996) in an adaptive mesh refinement method, to determine regions where 
refinement should be done. The other is the modification of the multiresolution method 
of Harten (1995) by converting it to a new, redundant, non-orthogonal wavelet. The 
wavelet sensor is then obtained by computing the estimated Lipschitz exponent of a cho- 
sen physical quantity (or vector) to be sensed on a chosen wavelet basis function. Both 
wavelet sensors can be viewed as dual purpose adaptive methods leading to dynamic 
numerical dissipation control and improved grid adaptation indicators. Consequently, 
they are useful not only for shock-turbulence computations but also for computational 
aeroacoustics and numerical combustion. In addition, these sensors are scheme inde- 
pendent and can be stand alone options for numerical algorithm other than the Yee et 
al. scheme. 


1 Introduction 

Efficient and accurate numerical simulation of fluid flows containing both sharp layers 
and turbulence are computationally very challenging. Resolving a wide range of length 
scales, including shock layers, and high shear mixings is time consuming and costly. 
Numerical methods of the total variation diminishing (TVD) type for shock capturing 
are too dissipative to be useful when turbulence is present. Higher order difference 
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department of Numerical Analysis and Computer Sciences, KTH, 100 44 Stockholm, Sweden. 
3 NASA Ames Research Center, Moffett Field, CA 94035 
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methods in conjunction with scalar or characteristics type nonlinear numerical dissipa- 
tions, without appropriate adaptive sensor control, behave similarly. Accurate methods 
such as spectral, spectral elements and high order spectral-like compact schemes for 
computing turbulence, break down when shocks are present. Although CPU inten- 
sive high order essentially non-oscillatory (ENO), weighted ENO (WENO) and discrete 
Galerkin schemes, generally exhibit less numerical dissipation than TVD schemes, nev- 
ertheless, their built in numerical dissipation still prevents the accurate capturing of fine 
scale turbulent structures without resorting to very fine grids. On the other hand, the 
switching mechanisms for multi-dimensional complex flow structures in hybrid schemes 
(e.g., switch between spectral element and ENO schemes) are extremely complicated 
and frequent activations of the ENO schemes are expected. As a remedy for this situa- 
tion, the artificial compression method (ACM) based filter scheme was proposed in 
Yee et al. [24], In this method one time step consists of one step with a fourth-order or 
higher accurate non-dissipative spatial base scheme. Often an entropy split form of the 
inviscid flux derivative (Yee et al. [25] and Gerritsen & Olsson [3] ) is used along with a 
post processing step, where regions of oscillation are detected and filtered by adding the 
numerical dissipation portion of a shock capturing scheme at these parts of the solution. 
The entropy splitting of the inviscid flux derivative is considered as a conditioned form 
of the governing equations. The idea of the scheme is to have the spatially fourth-order 
or higher non-dissipative scheme activated at all times and to add the full strength, 
efficient and accurate numerical dissipation only at the shock layers. Thus, it is nec- 
essary to have good detectors which flag the layers, and not the oscillatory turbulent 
parts of the flow field. It was shown in Yee et al. [24, 25, 26] that the ACM sensor, 
while minimizing the use of numerical dissipation away from discontinuities, consists 
of tuning parameters and is physical problem dependent. The objective of the present 
work is to develop adaptive numerical dissipation sensors that are improvements over 
the ACM sensor. 

Wavelets were originally developed for feature extraction in image processing and for 
data compression. It is well known that the regularity of a function can be determined 
from its wavelet coefficients [1, 13, 8] far better than from its Fourier coefficients. By 
computing wavelet coefficients (of a suitable set of wavelet basis functions), we obtain 
very precise information about the regularity of the function in question. This infor- 
mation is obtained just by analyzing a given grid function. No information about the 
particular problem which is solved is used. Thus, wavelet detectors are general, problem 
independent, and rest on a solid mathematical foundation. 

As of the 1990’s, wavelets have been a new class of basis functions that are finding use 
in analyzing and interpreting turbulence data from experiments. They also are used for 
analyzing the structure of turbulence from numerical data obtained from DNS or large 
eddy simulation (LES). See Farge [2] and her later work, and Perrier et al. [15]. There 
are several ways to introduce wavelets. One standard way is through the continuous 
wavelet transform and another is through multiresolution analysis, hereafter, referred to 
as wavelet based multiresolution analysis. Mallet and collaborators [8, 9, 10, 11, 12, 13] 
established important wavelet theory through multiresolution analysis. See references 
[20, 21] for an introduction to the concept of multiresolution analysis. Recently, wavelet 
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based multiresolution analysis has been used for grid adaptation (Gerritsen & Olsson 
[3] ) and to replace existing basis functions in constructing more accurate finite element 
methods. Here we utilize wavelet based multiresolution analysis to adaptively control 
the amount of numerical dissipation that is inherent in standard high-resolution shock- 
capturing schemes. With a proper choice of a set of the wavelet basis functions, we 
have a better control on the proper distribution of numerical dissipations leading to 
a more accurate simulation than the ACM sensor. The resulting wavelet sensors are 
readily available as more desirable grid adaptation indicators than those commonly 
used. It is well known that numerical dissipation is the main cause of wrong speed 
of propagation of discontinuous data in numerical combustion (LeVeque & Yee [7]) 
unless an order of magnitude adaptive grid refinement and time step reduction are 
used. Numerical dissipation is also a major stumbling block in efficient and accurate 
simulation of aeroacoustic problems. Consequently, the proposed wavelet based adaptive 
numerical dissipation control and grid adaptation indicator can be valuable to numerical 
combustions and computational aeroacoustics. In addition, this dual purpose adaptive 
method is scheme independent and can be a stand alone option for a variety of schemes 
other than what is discussed here. 

Section 2 reviews the Yee et al. [24] high order scheme employing ACM as a nu- 
merical dissipation sensor. Section 3 derives two multiresolution wavelet numerical 
dissipation sensors accompanied with scalar examples. Section 4 discusses the options 
in applying the wavelet sensor for highly coupled nonlinear systems of conservation 
laws. Numerical experiments for 2-D compressible Euler and Navier-Stokes equations 
are given in Section 5. 

2 High Order ACM Based Filter Scheme 


In vector notation the 2-D compressible time-dependent Euler equations in conservation 
form for a perfect gas can be written, in Cartesian coordinates, as 

Ut + F x + G y = 0, (2-1) 

where Ut = ^jf, F x = ^ and G y = ^ and the U, F, G, are vectors given by 

U= (p,pu,pv,e) T , 

F = ( pu , pir + p, puv , eu + pu) , G = (pv, puv , pv z + p, ev + pv) . 

The dependent variable U is the vector of conservative variables, and (p,u,v,p) T is the 
vector of primitive variables. Here p is the density, u and v are the velocity components, 
pu and pv are the x- and y-components of the momentum per unit volume, p is the 
pressure, e = p[e + ( u 2 + i> 2 )/2] is the total energy per unit volume, and e is the specific 
internal energy. For constant specific heats (calorically perfect gas) e = c v T , where c v 
is the specific heat at constant volume. 

The eigenvalues associated with the flux Jacobian matrices of F and G are (u, u, u±c ) 
and (v,v,v ± c), where c is the sound speed. The two u,u and v,v characteristics are 
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linearly degenerate. Hereafter, we refer to the fields associated with the u ± c and v ± c 
characteristics as the nonlinear fields and the fields associated with the u, u and v. v 
characteristics as the linear fields. 

Spatial Discretizations: The spatial discretizations of the ACM based filter scheme 
of Yee et al. [24] consist of two parts, namely, a base scheme and a filter. When 
numerical dissipations or filters are not used, the scheme consists of only the base 
scheme. If entropy splitting (Yee et al. [25] and Gerritsen & Olsson [3]) is used, the base 
scheme is applied to the split form of the inviscid flux derivatives. It is noted that the 
entropy splitting form of the inviscid flux derivatives can be viewed as a more conditioned 
form of the Euler equations for stability considerations. See Yee et al. [25, 26, 17] for 
details. Possible non-dissipative high order base schemes for F x and G y and the viscous 
terms (if present) are the standard fourth and sixth-order compact and non- 
compact central schemes for the interior grid points. 

The ACM Filter: For efficiency and ease of numerical boundary treatment, Yee et 
al. [24] proposed using filter operators whose grid stencils have a width similar to that 
of the base scheme. The filter operator consists of the product of a sensor and a 
nonlinear dissipation. Denote F,jj. as the discrete approximation of the inviscid flux 
F at (j Ax, A; Ay), where Ax and Ay are the grid spacings in the x- and y-directions and 
j and k are the corresponding spatial indices. Let the filter vector in the ^-direction be 
of the form 



= -R., i<5>* i. 
2 J+ 2 4+f 


(2.3) 


F* ! is the modified form of the nonlinear dissipation portion of the standard 
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numerical flux. For characteristic based methods, the quantity i? , i (with the k index 
suppressed) is the right eigenvector matrix of using Roe’s average state (Roe’s ap- 
proximate Riemann solver). G* k+1 is defined in the same manner. The elements of 

(with the k index suppressed), denoted by ((f3 + i)% are 



(2.4) 


4> l i is the dissipative portion of the high resolution scheme resulting from using a 

+ 2 ____ 

TYD, MUSCL with slope limiters, ENO or WENO scheme. Formulae for <f> 1 . i are 
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well known and can be found in the literature. See Yee et al. [24] for details and for 

a discussion of other possible filters. Here S l . , is the sensor and is a mechanism for 
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controlling excess numerical dissipation that is inherent in the dissipative portion of 
standard high-resolution shock-capturing schemes. 


For the numerical experiments to be presented later, we use the Harten and Yee 
upwind TVD numerical dissipation (Yee [22, 23]) 




)(9j + 1 + 9j) ~ ^{a l j+ i 


+ 7 


./• . Un / ’ '7 


(2.5) 
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t' + . = 


(sj+i -sj)/«J- + i «' + i^° 


S i + i = 0 


The . / = 1,2, 3, 4, axe the characteristic speeds of ! j£j evaluated at the Roe’s 
average state. The 5*. +i are elements of />’. 1 j. ( ^ Ty-;- 1 ./.• — Uj, k)- The corresponding 3*. +i , 
(b l . i and i? ■ , i using the MUSCL formulation are instead functions of the left and right 

states of f/. The function ip is an entropy correction to \a l 1 1. One possible form is 

3+ 2 

Harten & Hyman [6]. 

The ‘limiter’ function f/j used for the numerical experiment can be chosen as one of 
the following expressions 


= {<5'_i[( 5 hi )2 + ^ + 5 i+i [(S hi )2 + / 

= minmod(25' i, 2«' i , \(at\ , + a 1 ,)), 

J 2 J' 2 J' 2 J 2 

= 5 • max 0, min(2|5* 1 1, S ■ a l ._i ), min(|5* . . |, 

2 3 2 2 


( a j + i) +( a j_i) + 2 ^2 , 


25 • a 1 . 


Here 82 is a small dimensionless parameter to prevent division by zero and S = sign(5l , ) 
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The minmod function of a list of arguments is equal to the smallest number in absolute 
value if the list of arguments is of the same sign, or is equal to zero if any arguments 
are of opposite sign. 

The ACM Sensor: An artificial compression method (ACM) originally proposed by 
Harten [4] makes use of the gradient sensor to switch between a higher order scheme and 

a first order scheme. Yee et al. [24] use this gradient sensor as part of S l . 1 . In contrast 

3 ' 2 

to hybrid schemes that switch between spectral or spectral-like non-shock-capturing 
schemes and high order ENO schemes, the high order non-dissipative base scheme is 
always activated. For the ACM sensor, 5f +i = The parameter k is problem- 

dependent. The function 6 [. , is the Harten ACM gradient sensor. For a general 2m + 1 

3+2 

point base scheme, 


= max ( 6 


j— ra+l> *••5 


iyfi-iyji 

I a 1 . 1 1 + \a[. x I 
3+2 3 o 


Here the parameter p is an exponent > 1 and is not the “pressure p” in (2.1) and (2.2). 

For smooth flows in the absence of high shears, k can be very small. It is used to 
minimize spurious high frequency oscillation producing nonlinear instability associated 
with pure central schemes, especially for long time integration problems. Different 
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physical problems require different values of k because of the large variation in flow 
properties. The k value may vary from one characteristic wave to another, and from 
one region of the flow field to another region with different flow structure. Instead of 
varying k for the particular flow problem, one can vary p. For larger p. less numerical 
dissipation is added. Note that by varying p > 1 in the ACM sensor, one can essentially 
increase the order of accuracy of the filter. 

For a chosen grid spacing without grid adaptation, we would like to point out that 
neither the ACM sensor nor the wavelet sensor (to be discussed next) would be able to 
improve the accuracy at the shock and shear locations over the inherent shock-capturing 
capability of the nonlinear dissipation. The accuracy of the shock and shear is dictated 
by the chosen flux limiter of the nonlinear dissipation. The role of the sensors is to 
allow the full amount of numerical dissipation in shock and shear regions, and to limit 
the amount of numerical dissipation in regions immediately away from shock and shear 
locations and the rest of the flow field. Therefore, with a suitable sensor and flux limiter, 
one does not have to use CPU-intensive high order high-resolution shock-capturing 
numerical dissipation. The dissipation with sensor and flux limiter generally gives a 
slightly more accurate solution away from discontinuities but exhibits similar shock and 
shear resolution as second or third-order high-resolution numerical dissipation. 

Full Discretizations: If a multistage time discretization such as the Runge-Kutta 
method is desired, the high order non-dissipative spatial differencing base scheme is 
applied at every stage of the Runge-Kutta method. If viscous terms are present, they 
use the same order and type of base scheme as for the inviscid terms. There are two 
methods for applying the characteristic filter. Method 1 is to apply the filter at every 
stage of the Runge-Kutta step. Method 2 is to apply the filter at the end of the full 
Runge-Kutta step. For inviscid and strong shock interactions, method 1 might be more 
stable. 

If one desires a time discretization that belongs to the class of linear multistep 
methods (LMMs), e.g., trapezoidal rule or three-point backward differentiation, then 
the filter can be applied as a numerical dissipation vector in conjunction with the base 
scheme. The filter in this case is evaluated at U n for explicit LMMs. For implicit 
LMMs additional similar filters evaluated at the n + 1 time level might be involved. 
Alternatively, method 2 can be applied to LMMs as well. In this case, we apply the 
filter after the completion of the implicit time step. 

As an example, we illustrate the complete form of the schemes for Runge-Kutta 
methods with the filters applied at the completion of a full Runge-Kutta time step. 
Let U n+l be the solution after one full Runge-Kutta time step using a non-dissipative 
spatial base scheme. Then the solution at the next time level U n+1 is 


U "+ 1 = U ” +1 + ^ 

3 /\ 


F* i — F* i 

i+2 3 


+ 


At 




Here, the numerical filters F* i and G* , , i 

3 -^- 


A y [ 

are evaluated at U n+1 . 


( 2 . 10 ) 


6 



3 Wavelet Detection Algorithms 


The ACM sensor function is not entirely satisfactory since tuning parameters k and p 

are involved. We will develop here a wavelet based sensor function to replace k8 1 . ,, 

;+ 2 

which has the advantages of relying on a solid theoretical foundation, and minimizing 
the number of problem dependent parameters. Below we first describe the method in 
standard wavelet framework. Next we show how the Harten multiresolution method [5] 
can be used as a starting point for the derivation of a detection algorithm, on the same 
form. This modified Harten multiresolution description has the advantage of being more 
intuitive. 

There exists today a large body of results on determining regularity of a function 
from its wavelet representation. See for example [1], [12], [13]. The wavelet technique 
has been especially useful in simulation of turbulence, where it can be used as a data 
analyzing tool, extracting structures from the numerical solution, [2], [14, 15]. In [3] 
a wavelet based multiresolution analysis from [13] is used to determine points where 
mesh adaption should be done. We will, to some extent, follow the description in [3] 
but rather utilize the wavelet technique for numerical dissipation control. 

The wavelet decomposition of a function f(x) is obtained by taking the inner product 
with wavelet functions 'ip m , n {x). This leads to the wavelet coefficients, 

w m ,j = {/, = f f(x)ip mJ {x) dx, ^ 

m = . . . , —1,0, 1 . . . , j = 1,0, 1, . . . 

Here m is an index representing scale, and j is an index representing position. The 
set of basis functions is obtained by scalings of a single “mother wavelet” basis 

function 'ip(x) and is not to be confused with the ip 1 . i in (2.6). The construction of 

J~r 2 

ip{x) depends on the types of application. Among the many rich mathematical wavelet 
developments, one of the key elements of the present interest is the scaling of the mother 
wavelet. The other is the mathematical characterization of multiresolution scales by the 
so called Lipschitz exponents [1, 13, 12], 

There exist two major scalings leading to orthogonal and non-orthogonal wavelets. 
There also exist different scaling factors giving rise to different normalizations. The 
scaling we will use is 

ip m , j (x)=2~^((x-j)/2 m ). (3.2) 

With this scaling the resulting wavelets do not form an orthonormal basis, and this is 
sometimes referred to as a redundant wavelet decomposition. In addition, the mother 
wavelet ip{x) should have compact support, or rapidly decrease towards infinity, so that 
the wavelet coefficient w m j contains information about how much of scale 2 m is present 
in f(x) at the point x = j. Additional technical conditions on ip(x), for example that 
its integral is equal to zero, should be imposed. Under these conditions it is possible to 
compute f(x) from given coefficients w m . j . 

It is even possible to make the functions ip m j{x) an orthonormal basis for L 2 . In fact, 
more known results and applications are based on orthonormal wavelet basis functions. 
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In order to construct an orthonormal basis, the scaling of the wavelet function (3.2) is 
replaced by 

1 /; md (x) = 2- m ^(x/2 m -j) (3.3) 

so that the number of coefficients on the coarser scale is not as densely distributed. 
In practice, the scaling is restricted by the chosen grid size in numerical computations 
and consequently, cannot be very small. If rn corresponds to grid levels, this scaling 
(3.3) would restrict us to fewer coefficients on the coarser grids, which is not unnatural. 
For our case, however, since we will estimate the regularity of a function at all grid 
points using several levels of wavelet coefficients, we keep all level of coefficients at all 
grid points. We can then estimate the regularity of the function f(x) with the same 
accuracy at all grid points. With orthogonal wavelets we would have had a poorer 
estimate in between the coarse grid points. For a redundant wavelet decomposition, 
the scaling (3.2) is more natural and meets our requirement. Thus we will not use 
orthogonal wavelets. 

We are here mainly interested in computing the wavelet coefficients (3.1) and the 
information which we can obtain from them. In the wavelet based multiresolution 
analysis, there exist several theorems about the relation between the regularity of a 
function f(x) and its wavelet coefficients (3.1), see [1], [12], [2], [15], [8]. For example 
Theorem 9.2.2 in [1], which states that if ip is in C 1 and has compact support, and if 
the wavelet coefficients satisfy 


max | | < C2 ma (3.4) 

jeS(x 0 ,m,e) 

for all rn smaller than some limit, and for some e > 0, and some a, 0 < a < 1, then the 
function / has Lipschitz exponent a at x = Xo, 

\f{x)-f{x 0 )\<C\x-x 0 \ a (3.5) 

for all x in a neighborhood of xq . S(xo,m,e) is the extended domain of dependence of 
the wavelet function, 

S(xo,m, e) = {j : ip m j(x) ^ 0 for some x E (xo — e, Xo + e)} . (3.6) 

The assumptions on / and ip for similar results to hold vary between different refer- 
ences. Often it is required that / is in L 2 , bounded and continuous, and that ip is 
C 1 and compactly supported. Furthermore, the number of vanishing moments of ip(x) 
determines an upper limit on a. For example, if the first two moments of ip(x) vanish, 
the result is valid for 0 < a < 2. A kth vanishing moment means that f x k ip(x) dx = 0. 
The result in [13, 12], which is used in [3], instead gives the condition that ip{x) should 
be the ^'-derivative of a smooth function r](x) with the property 

i](x) > 0, / 7](x) dx = 1, lim r]^ k \x) = 0. (3.7) 

J rc^-ztoo 

Then the result is valid for 0 < a < k. A continuous function has a Lipschitz exponent 
a > 0. A bounded discontinuity (shock) is a = 0, and a Dirac function (local oscillation) 


8 



has a = — 1. Large values of k can be used in turbulent flow so that large vortices or 
vortex sheets can be detected. Although the theorem above does not hold for a negative, 
a useful upper bound on a can be obtained from the wavelet coefficient estimate. The 
work [8, 12] are good references for more detailed results on regularity estimates from 
wavelets. Before the discussion of the numerical computation of the wavelet coefficients 
we would like to emphasize that the Lipschitz exponent a in (3.4) and (3.5), which 
measures the regularity of a function f(x), holds the key to our wavelet based adaptive 
numerical dissipation control. 

3.1 Numerical Computation of the Wavelet Coefficients 

For practical computations, we cannot make the scale 2 m infinitely small. The smallest 
scale is given by the grid, which we normalize to m = 0. For the discrete case, the 
index m denotes the grid level, and the index j is the grid point index of the wavelet 
coefficients. The grid points are x = j, j = 1,2 ,N. For increasing m the support of 
'ipm.j ( x ) increases, which means that the scale becomes larger. For the present applica- 
tion, the regularity of the function that we want to analyze consists of numerical data 
obtained from numerically solving a system of nonlinear PDEs with a chosen numerical 
scheme (e.g. (2.10)). 

For a given grid function fj , j = 1,2,..., N that we want to analyze, we now 
describe how to compute numerically the coefficients w m j for m = 1 , 2 ,..., mo and j = 
1, 2, . . . , N. m o is the coarsest scale that we want to use. Of course, we could evaluate 
the coefficients directly from the definition (/, by numerical quadrature. This 

would be very expensive because the support of ip m j increases when rn increases, and 
the quadrature formulas would then involve sums where the number of terms approaches 
N on the coarsest scales. 

The basis of a computational scheme in obtaining the wavelet coefficients is the 
introduction of a so called scaling function (f>(x), which belongs to the class of smoothing 
functions r]{x) in (3.7) satisfying 

' 7 Q 

4>(x) = 2 dk4>( 2x — k) ip( x ) = 2 Ck(f>( 2x — k). (3.8) 

k=—p k=—p 

Here, the scaling function (f>(x) has no association with the scaling of the mother wavelet 
ip{x) by 2 m . Also, the scaling functions <p(x) and <f> m j(x) are not to be confused with the 

nonlinear numerical dissipation (f>[. i in (2.4). Only a few dp. and cp should be nonzero, 

3 '2 

in order to obtain an efficient computational method. 

The functional equations (3.8) is closely associated with the multiresolution anal- 
ysis [20, 21, 8]. In the multiresolution analysis, a sequence of closed subspaces Vj of 
L 2 (M), j E Z of different resolution is defined, 

. . . C V 2 C Hi C F 0 C V-x C V-2 ■ ■ ■ 

with the important property that if a function h(x) E V rn . then h(2 m x) E Vo, so that 
all spaces are scaled versions of Vq. Here M and Z are the spaces of reals and integers 
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respectively. Furthermore, we require that 


|J V m = L 2 (W) and p) V m = {0}. 

m£ Z 

In other words, \J me v V m is dense in L 2 ( R). The spaces V m can be thought of as 
representing all scales down to the mth scale. We also require that Vo is invariant under 
integer translations, so that a basis on the form {(/){x — j ) } j can be found. With the 
definition 0 rn . ? = 2~ m / 2 0(a:/2 m — j), we then obtain {4>m,j}j as a basis for V m . With 
these definitions the scaling relation (3.8) for 0 just means that we expand the function, 
e.g., 0o,o in the basis 0_i,j, which can be done since Vo C V-\. 

Let W m be the orthogonal complement of V m in V m j , so that 

Vm - 1 = W m V mi 

where the symbol 0 stands for direct sum. It can then be shown that a basis {ip m ,j{x)} 
for W m can be found, where i /> m j = 2~ m l 2 'ip{x /2 m — j). The spaces W m are orthogonal, 
since W m C V rn - 1 , and V' m _ i TII-y a _ i . We can write L 2 as a direct sum of W m spaces, 

L 2 (R) = 0 W m 

m 

so that the ip m j form an orthogonal basis for L 2 . This is an orthogonal wavelet basis 
for L 2 . The scaling relation for -0, can then be thought of as an expansion of 0o,o i n the 
basis 0-i, j, which can be done since Wo C V_i. 

If a scaling function <f>(x) satisfying (3.8) can be found, we obtain a fast method 
for the evaluation of the wavelet coefficients, independent of whether the formula was 
derived from an orthogonal wavelet basis, or was found in any other way. The compu- 
tational scheme becomes 


(/> 0m, j) — dk (/j 4 > m—l,j+2 m ~ 1 k ) 

(/ 5 0m, j} X^A-=— p (/i 0m— %,j+2* a ~ 1 k) 


where we have defined 0 mJ (x) = 2 m (p((x — j)/2 m ). To derive (3.9) from (3.8), insert 
the argument (x — j)/2 m in the place of x and multiply by 2~ m to obtain 


k=—q 


X -J 


~k) = 


1 p 

— T 


k=—q 


x ~ j ~ 2 

2'm-l 


m— 1 } 


from which we can identify the relation 


p 

4 > m,i(x) = 'y ^ d^(p m _i j + 2 m - 1 k{x)- 
k=—q 


Multiplying this equation by f(x) and integrating, gives (3.9) for 0, the relation for 
0 follows from similar computations. If the coefficients on the finest grid (/, 0o,j) are 
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given, we can compute {/, 'ipij), (/, by using the formulae (3.9). Repeating the 
procedure gives for successively increasing m. The formulae (3.9) are a pair 

of difference operators acting on a grid function fj. We can express the algorithm as 
follows. Introduce the grid operators 


Afi 

Dfj 

and its '/nth level expanded versions 


— EL- P dkfj+k 

p ^kf j+k 


Amfj — E k=—p dk f i+2 m k 
Dm.fj = Efc— — p c kf j+2 m ki 


(3.10) 


(3.11) 


where the integers p and q are related to the chosen ip(x) and 4>{x). The computation 
(3.9) can then be written as 


/<•"> = </,.w = -wj”- 11 (312) 

W m j = ■ 

Starting from the grid function /j 0 ^ = oj), (3.12) is applied successively to obtain 
grid functions for rn = 1, 2, 3, ... , itiq. The initial /j 0 ^ is found by numerical quadrature. 

The computation /j m - > = A, m |/j m 1 ^ for a three point operator A is outlined in Fig. 3.1. 
These sequences are sometimes referred to as the impulse response of low pass and 
high pass filters. The computation w m ..j = D rn i/j™ 1 ^ follows a similar pattern, but 
possibly with a different stencil width. 



The mth level of wavelet coefficients can be written as 


— {f i — Dm—lAm— 2-^-m — 3 • • • ^0 fjt ITl — 1,2,.... (3.13) 

Once the coefficients and Cj. are determined the computation is a very standard 
application of grid operators. In practice, we only use mo = 3 to 5. To be able to 
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compute up to the boundary, we use one sided versions of the given operators. This 
seems to work well in practice, although it is not covered by the wavelet framework 
described above. Study on the effect of wavelets on an interval using the appropriate 
boundary wavelets [8] for our application is a subject of ongoing research. 

To determine the initial coefficients, {/, fo,j) from fj , we make a numerical approxi- 
mation of the integral. If the support of f(x) is small, we can set {/, <fioj) ~ fj f 4>{x) dx. 

3.2 Approximation of the Lipschitz Exponent 

After we compute the wavelet coefficients, the next step is to compute the Lipschitz 
exponent a, the key for our present development. Unfortunately, it is not possible to 
obtain a exactly. There exist in the literature different methods for approximating a. 
The method of approximating a in [3] is based on the theorems in [13], which involves 
tracing a maximum curve among the wavelet coefficients to points of singular behavior. 
We choose here, instead, to base our method on Theorem 9.2.2 of [1] as described 
previously. This means that we do not trace maximum lines. Instead the procedure 
below is applied at all grid points. Note that the results in both [1] and [13] are valid 
for functions of a continuous variable, so there is some freedom of interpretation when 
applying them to functions defined only at grid points. 

After we have computed the wavelet coefficients, we first form the maximum over 
the domain of dependence, 


r m,j — , I (fitpmj+k) I (3-14) 

k=—2 m p,2 m q 

where the nonzero coefficients are enumerated from —p to q. We estimate the Lips- 
chitz exponent by a least squares fit of a line to the equation 

log 2 r m ,j = am + c. (3.15) 

The slope gives an estimate of a at the point x = j. A discontinuity is characterized 
by a = 0. Standard numerical centered difference approximations have problems when 
a is small. Usually existence of several derivatives is required for high order difference 
formulas to be accurate. 


3.3 Detectors from the B-Spline Wavelet Basis Function 


Developing the best suited wavelets that can characterize all of the flow features might 
involve the switching or blending of more than one mother wavelet ip{x) and scaling 
function (f>(x). The mother wavelet function used in [3] and described in detail in [13] 
meets some of our requirements. It is obtained from second order B-splines. 


<ip(x) 


0 

— 2{x — l) 2 
—Ax(A — x) + 2x 2 
— Ax(l + x) — 2x 2 
2(x + l) 2 
k 0 


x > 1 

1/2 < x < 1 
0 < x < 1/2 
-1/2 < x < 0 
-l<x< -1/2 
x < — 1 


(3.16) 
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For this wavelet (3.16), there exists a scaling function, given by 

1 0 x > 2 

— 2) 2 1 < x < 2 

—(x — 1/2) 2 +3/4 0 < a; < 1 . (3.17) 

\{x + l) 2 — 1 < a; < 0 

0 a; < — 1 

If we apply the Fourier transform, h(£) = f e~ lx ^h{x) dx, to the relation 

4>{x) = 2 dk<t>( 2x — k ) 
k 

we obtain 

m= bottom/ 2) (3.18) 

where &o(£) = Ylk^k elk ^- From a given trigonometric polynomial &o(£)> we can fi nc l 
scaling functions by iterating (3.18), 


m= 


m = 1 

or conversely, for a given scaling function, we can try to find the polynomial &o(£) from 
(3.18). Taking the N degree B-spline basis function as the scaling function gives 

where s is 1 for N even, and 0 for N odd. From this it is not hard to verify that 

b 0 {O = e 2 c..s v -'^/2 


and we get the coefficients d^ from this trigonometric polynomial. 

The normalization is such that the integral of the scaling function above is equal to 
one. The functions above are standard, and can be found in [1], The scaling function 
differs by a shift from the scaling function used in [3], but the important relations 


4>{x) = \<j>{ 2x + 1) + i<j>(2x) + \(j>( 2x - 1) + \(j){ 2x - 2 ) 
4'(x) = 4>( 2x + 1) — <^(2a:) 

hold, and give the grid operators 

Afj =(fj-i + 3/i + 3/ j+1 + / i+2 )/8, j = 2, N — 2 
Dfj (/, i fj)/2 J 2 V. 
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Scaling Function 



Wavelet Function 



Fig. 3.2. Scaling function. 


Fig. 3.3. Wavelet function. 


Here p = 1, q = 2, d- \ = 1/8, do = 3/8, d\ = 3/8, <ft = 1/8, c_i = 1/2, co = 0, 
ci = 1/2, and C 2 = 0 for (3.11). For example, for mo = 3, the computation of /j m ^ 
involves a grid stencil of 13 points. Note that the wavelet stencil is not symmetric. In 
general, formula (3.9) shows that points from p2 m °~ 1 to —q 2 mo_1 are involved in the 
computation, giving a stencil of totally (p + q) 2 m ° -1 + 1 points. The scaling function 
4>{x) and the wavelet function ip(x) for the B-spline wavelet are shown in Figs. 3.2 and 
3.3 respectively. 

For boundary operators we use 


Dfi = (ft- 3/ 2 + 2ft)/2 
Ah = (7ft — 3/ 2 + 5ft — / 4 )/8 
Afw-i = (7ftv — 3/at-i + 5/a^_2 — /at-3)/8 
Afw = (25/ \ - 37/ v 1 + 27/jV- 2 - 7/. v 3 )/8. 

The coefficients are determined so that the action on smooth functions is close to the 
interior operators. The simpler formulae 


Aft = (ft — ft)/2 
Aft = (7ft — 3ft + 5ft — ft)/8 
AfN-l = (/at + 3/aT-I + 3/aT-2 + /at-3)/8 
AfN = (7/at — 3/at-i + 5/at-2 — /a f -3)/8 

have also turned out to work well in practice. We will later show some experiments 
with these operators. The computation consists of the following steps. Use (3.20) with 
(3.21) or (3.22) as the operators (3.10). Compute from (3.12), and estimate the 
regularity by (3.14) and (3.15). 
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3.4 Detectors from Converting Harten’s Multiresolution Scheme to 
Redundant Wavelets 

This section describes a procedure to modify Harten’s multiresolution method ([5]) 
into redundant non-orthogonal wavelets and derive an alternate wavelet detector which 
in some sense is slightly simpler than the B-spline wavelet sensor. The multiresolution 
method by Harten is a way to speed up computations by using a wavelet decomposition. 
See Sjogreen [19] for a study of Harten’s multiresolution scheme. We next give a brief 
description of the method. 

Consider the approximation of a PDE in one space dimension, on the uniform grid 
Xj = j Ax, j = 0, 1, . . . , N. The numerical solution is given by (/o, / 1 , • • ■ , /a-), with fj 
an approximation of the solution at xj. Introduce L levels of successively coarser grids, 


G k (xq } X 2 , . . . , ) k 0 

Let Xj denote grid point j on grid G k . Then Xj = j Ax, and x k = Xj 2k = j 2 k Ax. Let 
N k denote the number of points in G k ■ Then N k = 2 L ~ k NL. We let f k denote the 
numerical solution of a PDE at the point x k . 

Assume that the solution is given on grid G k , and that we want to find it on the 
finer grid G k _ 1 . For the even numbered grid points we have 

/*; : ,/;f- j = o,i,...,N k . 


To find the solution at the odd grid points, we let I(x,f k ) interpolate f k on G k , such 
that I(xj,f k ) = f k . We then have the approximation f k f-\ = f k ) of 

The interpolation error is d k = ~ f‘i j-\ ■ Thus with knowledge of f k and d k we can 

reconstruct the solution on G k ~ 1 . We call ( f k ,d k ) the multiresolution representation 
of f k ~ 1 . Note that the vectors f k ,d k together contain the same number of elements as 
does In summary, we switch between the representation /* -1 , and ( d k . f k ) by the 

forward transformation 


fk fk — 1 

Jj J2j 


:= / 


k - 1 

2j-l 


I(x 


k - 1 
2j-V 


f k ) 


(3.23) 


and the backward transformation 


fk — 1 . fk 

J 2 j —Jj 

d) ■ I[.v^\.f k ). 

This is inexpensive if I(x k ~\, f k ) is a straightforward linear interpolation operator. 

We transform consecutively on all grids 

f° -> (rf 1 ,/ 1 ) (d\d 2 ,/ 2 ) ...(d\d\...,d L J L ). 

The vectors (d 1 , d 2 , . . . , d L , f L ) are the multiresolution representation of /°. The inter- 
polation errors d k contain information about the smoothness of the solution. If d k is 
large at some grid points, this indicates that the solution is non-smooth there. 
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The multiresolution method does not increase the number of coefficients, and is 
therefore not redundant. We here propose instead to use a redundant form, obtained 
from computing coefficients d/j at all points on the fine grid. To do this, we treat the 
fine grid as consisting of two coarse grids. The even points are interpolated from the 
odd points, and vice versa. The interpolation errors are used for wavelet coefficients. 
The following formula then replace (3.23), 

f2j-n4j l J2r- 1) 
ft,-, ■■= n4j-v flf') 

By expressing the wavelet coefficients as interpolation errors, it becomes intuitively clear 
why they should be small when the function is regular. 

We express the computation of the level k coefficients from the level k — 1 coefficients 
in operator form as 

f} = A fr 

d k = fk-i _ Af k-i 

where Afj is an averaging operator, coming from the interpolation formula. We assume 
that the same I(x, f) is used for the odd and the even points. The computation is done 
for all j on the fine grid. At the boundary points, we use one sided versions of the 
operator A. Note that (3.23) is exactly the same form as (3.12), if we define D = I — A, 
with dj as w m j, and f 1 - as /j m ' 1 in (3.12). 

There is more than one choice for the interpolation function. See Sjogreen (1995) 
for a discussion. The exact form of the method for the computations in this article will 
be 

Afj =(/,-i+/,+ 1)/2 J 2 V I 

Dfj = fj — Afj j 2 V I. 

The above choice was made in order to have a simple and efficient method. The stencil 
is narrower than for the B-spline formulas that were given previously. With the formula 
above, we also get a symmetric stencil, which is more natural if the other parts of 
the computation, such as difference approximations of PDEs are done by symmetric 
formulas. Furthermore, symmetry makes periodic boundary conditions somewhat easier 
to implement. Note that the absence of symmetry for either the scaling function or the 
wavelet can lead to phase distortion. This can be shown to be important to signal 
processing applications. 

We now have a method in the same form as the wavelet detector from the previous 
subsection, but derived in an intuitive way and using a more narrow grid stencil. For 
example, for mo = 3, the computation of /j mi involves a grid stencil of 9 instead of 
13. Using the above formula, we proceed with the Lipschitz exponent computations as 
before. 


16 



3.5 Two Space Dimensions 

Wavelet functions can be defined in different ways in two dimensions. One way is to 
use the product basis functions, 

'4 > m,n,j,k(x) = tfm.,jix)'fn,kix) 

giving two scales rn and n at each point (j. k ) at the plane. It is possible to instead use 

***(*) = - i)/2 m . (v - *)/ 2”) 

where only one scale is present at each grid point. ip(x,y) is a two dimensional version 
of ip(x). 

The computation of multi-dimensional wavelets is quite expensive, especially in 3- 
D. A simple minded efficient way is to evaluate the wavelet coefficients dimension by 
dimension. This means that we get two set of wavelet coefficients wt, , , and vff 

''*’ 5 . 75 '*' ''* 5.7 5 ™ 

where now (j,k) is the position and rn is the scale. The precise definition is 

W m,j,k = J fix, k)lp md {x)dx 
W m,j,k = I fUMmMdy. 

For simplicity in notation, here f{x,k) means the function “/” to be sensed in the re- 
direction with a fixed k grid index in y. Later, our numerical method will have some 
terms evaluated as finite differences in the re-direction and some which are evaluated 
in the y-direction. We then use the wf n ■ k coefficients for the re-direction computation, 
and the y-coefficients for the y-direction computation. 


3.6 Numerical Experiments with 1-D Test Cases 


We here test the detector algorithms on one dimensional examples, where the regularity 
of the functions is known. We will compare the redundant wavelet of Harten’s mul- 
tiresolution method (3.24), referred to as RH- wavelet, and the B-spline wavelet (3.20) 
referred to as BS-wavelet. A typical example from the wavelet literature is the function 
given in Fig. 3.3. This function has a steep layer, a jump (« = 0), a sharp peak which 
on the grids used is underresolved and could be understood as a Dirac pulse. A smooth 
bump is also present. It has the form 


fix) 


( 0 


< 


l/( e - 15 (x+ 3 /4) + j) 

1 

0 


1 -20|®- 1.95 
0 

e -15(x-7/2 f 


y o 


x < — 1 

-1 < X < -1/2 
-1/2 < rr < 1/2 
1/2 < x < 1.9 
1.9 < x < 2 

2 < x < 3 

3 < x < 4 
x > 4 


(3.25) 


A good detector should flag the discontinuity, and the possible spike. 
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Figure 3.4 shows results from using the ACM sensor k6 1 . , in (2.8) and (2.9) with 

3 ' 2 

p = 1 on the function in Fig. 3.3. We here define 5*. + j as f j+ \ — f j. The cases k = 1 
and k = 0.5 are shown. Note that k is more relevant if it is a function of the dependent 
variables, and/or when more than one wave is involved as in systems of hyperbolic PDEs. 
Here k is just a scaling parameter. The value 1 (0.5) corresponds to non smoothness, 
and the value 0 corresponds to a smooth linear solution. The function was evaluated at 
300 grid points. We observe that the detector is unable to distinguish a corner (a = 1 
for the wavelet detectors) from a discontinuity (a = 0 for the wavelet detectors) since 
the ACM sensor is a single scale detector. Furthermore, the smooth maximum also 
triggers the detector. Figure 3.4 indicates that numerical dissipation with the indicated 
amount will be used for all of these nonzero regions. 



Fig. 3.3. Function used for testing. 



Fig. 3.4. ACM sensor, k = 1 (left) and k = 0.5 (right). 


Next, the BS-wavelet operators (3.20) are used to compute four levels of wavelet 
coefficients for the function given in Fig. 3.3, evaluated at 300 grid points. A least 
squares fit to the line (3.15) is done. The resulting local regularity exponent a is 
plotted in Fig. 3.5. If some wavelet coefficients are equal to zero, a is set equal to one, 
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since this corresponds to maximu m regularity and one is the highest a this particular 
wavelet is able to predict. The initial coefficients o,«) were evaluated by numerical 
quadrature from function values at the grid points. In this computation a is equal to 
0.13 at the jump, and —0.46 at the spike, and is close to or above one for the rest of 
the function. 

Figure 3.6 shows the same computation, but with the RH- wavelet decomposition 
using operators (3.24) instead. The same function as above is decomposed into four 
levels of coefficients, and a least squares fit is done to find a. The predictions are of 
similar type as for the BS-wavelet, but give a somewhat lower value of a. The RH- 
wavelet gives the a value of 0.0 for the discontinuity, and -0.65 for the spike. Here we 
made a numerical quadrature formula for the initial coefficients, (/, 4> o,n) = (fj+fj+ 1 )/ 2 . 
This was necessary to do in order to have an exactly correct result for step functions. 
In [13], [3] scaling factors, which change with the scale level, are used on the wavelet 
coefficients in order to renormalize for step functions. 

Readers are cautioned not to compare Fig. 3.4 with Figs. 3.5 and 3.6 in a straight 
sense. The ACM sensor k6 1 i has no one-to-one correspondence with the Lipschitz 

exponent a, even in the actual implementation of the wavelet sensor in the numerical 
scheme. Figure 3.4 shows the amount of numerical dissipation needed for the Yee et al. 
scheme to capture a solution like (3.25), whereas Figs. 3.5 and 3.6 show the regularity 
of (3.25). The amount of numerical dissipation determined by the wavelet sensors does 
not come into play until later. In a loose sense, Figs. 3.4 - 3.6 illustrate the fact 
that the ACM sensor is a single scale detector as opposed to the chosen multiple scale 
multiresolution wavelet sensors which are capable of detecting all of the four features of 
the function (3.25). 



Fig. 3.5. BS- wavelet a estimate. 



Fig. 3.6. RH-wavelet a estimate. 


As a second example, we investigate the capability of the detectors to predict the 
exponent of the function \x\^, where 0 < [3 < 1. The Holder (Lipschitz) exponent a 
obtained from the BS-wavelet decomposition (3.20) with four wavelet levels is shown in 
Fig. 3.7 as a function of the exact exponent (3 where 300 grid points were used. The 
same quantity computed by the RH-wavelet decomposition (3.24) is shown in Fig. 3.8. 
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In the figures the exact exponent is also plotted as a line with slope one. We note that 
the prediction is fairly good in the interval 0.5 < (3 < 1, but that the computations tend 
to overestimate the regularity at the lower end of the interval. This is probably due to 
poor resolution of the singularity. 




3.7 Experiments with orthogonal wavelets 

To contrast the behavior of the non-orthogonal wavelets, we perform numerical experi- 
ments with an orthogonal wavelet. We derive a computational procedure for computing 
orthogonal wavelet coefficients from a relation 

(f>(x) = \[2 h/ c (f)(2x — k) (3.26) 

k=—p 

and make the scaling (f> m j(x) = 2~ m l' 2 (j>(2~ m x — j) . It is possible to show that a similar 
formula holds for the wavelet function, 

<i 

ip{x) = \P2 ^2 9kH 2x - k), 

k=—p 

with fip = (— 1 ) k h_p + i. The wavelet basis functions are then defined as ip m j(x ) = 
2- m/ V(2- m x — j). The formulae for the coefficients (/, 4>m,j) and (/, 'ipm.j ) now become 

(/ ) $m,j) = YH-- p dk (/i (pm-l,j+k) (3 27) 

(fkPm,j) = 'l2 < l=-p C k (f i ( Pm—l,j+k) 

instead of (3.9). These coefficients can be thought of as existing on a dyadic grid, with 
grid points 2 m j on scale rn. The grid for orthogonal wavelet computation is outlined in 
Fig. 3.9. Examples of orthogonal wavelets are given in ([1]). We have here chosen to 
implement the so called Daubechies wavelets, D(N ), ( see p. 194 in [1] ). These wavelets 
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have compact support and are orthogonal. The TVth wavelet of the family has a support 
of length 2 N — 1. The computations shown below have been done with N = 3. 


m = 2 


m = 1 


Fig. 3.9. Dyadic grid 


m = 0 
3 


Table 3.1. Coefficients for orthogonal wavelet D( 3). 


k 

hk 

0 

0.3326705529500825 

1 

0.8068915093110924 

2 

0.4598775021184914 

3 

-0.1350110200102546 

4 

-0.0854412738820267 

5 

0.0352262918857095 


Scaling Function, Daubechies wavelet 



Fig. 3.10. Scaling function. 


Wavelet Function, Daubechies wavelet 

2 
1.5 
1 

0.5 
0 

- 0.5 
-1 

- 1 . 5 ' * * * * * 1 

- 3 - 2-10123 

x 

Fig. 3.11. Wavelet function. 



The coefficients h/. for this wavelet is given in Table 3.1, where the limits in the 
summation in (3.26) are p = 0 and q = 5. The scaling function and the wavelet function 
are shown in Fig. 3.10 and Fig. 3.11 respectively. 

We can now compute just as before, but on the dyadic grid in Fig. 3.9, instead of a 
standard grid. We compute the wavelet coefficients by the formula (3.27). To estimate 
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the regularity at a point j, we form the maximum over the wavelets coefficients, 

r m ,j = max | (/,Vw) | 
k 

where k ranges over the values for which the support of 'ip m ,k{x) contains j. Here 
fewer coefficients than in the non-orthogonal case are, in general, used since they are 
distributed more coarsely on the dyadic grid. Once the quantities r m j are found we 
proceed as previously, remembering that the different scaling in the orthogonal case 
makes the wavelet coefficients decrease as 2 (a + l / 2 ) m with scale for a function with 
regularity exponent of a. In the non-orthogonal case the behavior was 2“ m . Figure 
3.12 shows the result from the above procedure applied on the testing function (3.25). 



X 


Fig. 3.12. a estimate from orthogonal wavelet D (3). 

While the orthogonal wavelet is able to detect all four features of /, it gives a entirely 
different a values, especially, the shock and the Dirac-like pulse. More experiments with 
estimation of regularity can be found on pp. 301-304 in [1], 

4 Wavelet Detectors in the High Order ACM Filter Scheme 

In the previous section we described how to estimate the regularity (Lipschitz) exponent 
a numerically. However, we have not yet discussed how to use the information contained 
in the Lipschitz exponent. One possibility would be to let the order of accuracy of the 
numerical scheme adapt to the regularity of the function. Other possibilities are either 
to use an approximation adapted to the regularity of the solution similar to the (h,p) 
finite element methods, or to integrate the wavelet sensor into a difference scheme, by 
using the Lipschitz exponent a in the same way as the switching quantity used in the 
high resolution shock-capturing schemes. These are subjects of ongoing investigations. 
In this paper, we discuss how to use the Lipschitz exponent to dynamically control the 
numerical dissipation. We concentrate on the simplest method of improving the Yee 
et al. ACM based filter scheme. The straightforward procedure is to insert the new 
estimator into the ACM based filter scheme to replace the ACM sensor. This boils down 

to how to adaptively switch on or gradually transition to the numerical dissipation (f>[. i 

J '2 

in (2.4). 
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For the numerical experiments presented in Section 5, the wavelet sensor is obtained 
by computing a vector of the approximated Lipschitz exponent of a chosen vector 
function to be sensed with a suitable multiresolution non-orthogonal wavelet basis 
function. Here, “vectors or variables to be sensed” means the represented vectors or 
variables that are suitable for the extraction of the desired flow physics. The variables to 
be sensed can be the density and/or pressure, the characteristic variables, the jumps in 

the characteristic variables a 1 it or the entropy variable vector W (Gerritsen & Olsson 

3' 2 

[3], Yee et al. [25]). The choice on the type of wavelet basis functions and their scaling 
functions depends on the types of feature that we want to extract or detect. For the 
test problems to be presented, we would like the wavelet basis function and its scaling 
function to be capable of detecting shocks, shears, spurious oscillations and turbulence. 
The BS and RH wavelets are two possible choices. Study on the optimal choice of 
wavelet basis functions and their scaling functions is a subject of ongoing research. 

For the wavelet sensor, the sensor <Sj +1 / 2 in (2.4) can be defined as 

Sj~ I 2 = T ( a j' + 1/2) (4-1) 

where «j +1/ / 2 ^ ie estimated Lipschitz exponent of the Ith characteristic component 
with I = 1,2, 3, 4, the four characteristic waves. r(a) is a sensing function which de- 
creases from r(0) = 1 to r( 1) = 0. Note that the Ith component of the estimated 
Lipschitz exponent is not to be confused with the jump in the Ith characteristic 

variables a 1 . , in Section 2. 

■?+2 

If we instead base the exponent estimate on point centered quantities, we will use 
the sensor function 

S j+i /2 = max(r(«‘),T(a' +1 )) (4.2) 

and if the exponent estimate is based on other quantities than the characteristic, e.g., 
density and pressure, we use the switch 

%i/2 = maxSj+i/2 ( 4 -3) 


where the maximum is taken over all components of the waves used in the estimate, 
and which thus is the same for all characteristic fields. 

The function r(a) should be such that r(0) = 1, and r( 1) = 0 or a smooth transition 
between 1 and 0. Three options are considered. 


, , [ 1 a < an 

T(a) = < 

[0 a > ao 

r(a) = 5 + 4 arctan if («o 
r(a) = max{0, min[l, (a — 


— a) 

l)/(a 0 -l)]}. 


(4.4) 


Here, «o is a cut off exponent to be chosen, and if is a constant, which we have tried 
with values in the interval [200, 500]. Alternatively, one can integrate the actual a value 
into the sensor function instead of using the same amount of numerical dissipation at 
the cut off exponent, especially to be used as grid adaption sensor as well. 
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After some experimentation, we have found that switching on the dissipation at the 
grid points where a < 0.5 works well, i.e. , taking 


1 a < 0.5 

0 a > 0.5 


(4.5) 


In fact the method does not seem to be very sensitive to the exact value of cut off 
«o, (for 0.4 < «o < 0.6) for all the test cases considered. Furthermore, the same cut 
off value, 0.5, works well for all problems we have tried in Section 5 (except for the 
vortex convection case, where «o = 0.0 to be used in conjunction with entropy splitting 
[25]). Experiments with smoothed step functions do not give very different results. To 
distinguish the high order ACM based filter scheme, we referred the scheme discussed 
in Section 2 using the wavelet sensor (4.1)-(4.2) as the wavelet based filter scheme. 

We would like to point out that the simple minded wavelet sensor (4.4) or (4.5) does 
not make full use of the Lipschitz exponent. Regardless of the amount of information 
on the Lipschitz exponent used to design the wavelet sensor, the sensor is not really 
parameter free. Unlike the ACM sensor, the parameter involved in the wavelet sensor 
is, however, not arbitrary. There are theorems and guidelines on what values of the 
Lipschitz exponent to be expected for various features of the function to be sensed. 


5 Numerical Experiments for 2-D Compressible Euler and 
Navier-Stokes Equations 


To illustrate the performance of the wavelet sensor, the same three perfect gas test cases 
with distinct flow properties as in Yee et al. [24, 25] are used. The first is inviscid and 
the last two are compressible full Navier-Stokes computations. The three test cases are: 
(1) a 2-D inviscid horizontally convecting vortex with periodic boundary conditions 
(BCs), (2) a 2-D vortex pairing in a time-developing mixing layer with shock waves 
formed around the vortices, and (3) a 2-D shock wave impinging on a spatially evolving 
mixing layer where the evolving vortices must pass through a shock wave, which in turn 
is deformed by the vortex passage. Figures 5. 1-5.3 show the schematic, flow conditions 
and the computational domains of the three test cases. 
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Isentropic Vortex Evolution 

(Horizontally Converting Vortex, vortex strength 0=5) 


Freestream : 

(ttoojWoo) = (1)0); Poo ~ poo 1 

JC: Perturbations are added to the freestream (not in entropy) 

(Su,Sv) = — e^(-(y - yo),(z - *„)) 
#r = ~ 1 )^ 2 e i-r a 

877T 2 

i 2 = (x-x o y + (y-yof 

Computational Domain & Grid Size: 


0 < * < 1 0 & —5 < y < 5 
80 x 79 Uniform grid 
Periodic BC in * & y 


Initial Vortex, Density Contours 



(*»« i y»o ) 
center of vortex 


Fig. 5.1. Test problem 1, isentropic vortex convection. 


Vortex Pairing in a Time-Developing Mixing Layer 

(M c =0.8, Re=1000, T R =300K, Prandt! #=0.72) 


Normalized with vorticity thickness: 


6 U = 


«1 -«2 


(du/dy) mam 

T & c: determined by assuming constant stagnation enthalpy 

1C: 

u-| a 1> u 2 =-1, T-j =T 2 , 


Inital shear profile: U = 0.5 tanh (2y) 

Crocco-Busemann: 
c 


■ = «; + Vt-f — *> 



Initial perturbations: 
2 


»' = cos ( 2lr ^ x / L * + L x = 30,6 = 10 

02 = 0.05, = — x/2 

BC : Periodic in x, slip walls in y = 0-01, <h = — w/2 

Grids: 


L y smh(b y T]) 

»=f-S5(M' = 100 ' = 3 ’ 4 

uniform in * 


Fig. 5.2. Test problem 2, vortex pairing. 
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Shock Impingement on a Spatially-Developing Mixing Layer 

(M c =0.6, Re=500, TR=300K, Prandtl #=0.72, oblique shock angle = 12°) 


Convective Mach #: , , «i - 

Me — 


1C: 


Cl + c 2 

Ui=3, u 2 =1 , equal p & stagnation enthalpies for the 2 streams 

Initial profile: u = 2.5 + 0.5 tanh (2y) 

Inflow perturbations: 

2 

.' = X> cos(27t kt/T + 4>u) exp(-i/ 2 / 

h=l 

period T — 30/2.68, b = 10 
ai = 0.05, <f>i = 0 
a 2 = 0.05, <f >2 = 7r/2 

B£: 

Oblique shock wave forced by upper BC (impacts the shear layer at x=90) 

supersonic outflow 

lower boundary specified as slip wall 

Grids: 



L y sinh(r/) 
2 sinh(l)’ 
uniform in 2 ! 


y = 


L * = 200, L y = 40 


Fig. 5.3. Test problem, 3, shock impingement on a mixing layer. 

In order to compare results with Yee et al. [24, 25], we use the same time and spatial 
base scheme as Yee et al. The classical fourth-order Runge-Kutta time discretization, 
and the non-compact central spatial schemes with the same order of accuracy and type 
of base scheme for the inviscid and viscous terms (if viscosities are present) are employed. 
The filters are applied at the end of the full Runge-Kutta time step. Roe’s [16] average 
states are used in (2.3), along with the Harten and Yee ([22, 23]) second-order upwind 

TYD dissipation portion (2.5)-(2.7) for (p[. i in (2.4). From numerical experiment, 

3 ' 2 

limiter (2.7a) produces the best result for test case 1 and limiter (2.7c) for test cases 
2 and 3. Results shown in Section 5.2 reflect these choices. For the ACM sensor the 
parameters p and rn in (2.8) and (2.9) are set to 1 and a small value of 10~ 6 is added to 
the denominator of (2.9) to avoid an extra logical statement for the ACM sensor. For 
the wavelet sensor, the cut off Lipschitz exponent «o is 0.5 for test cases 2 and 3 and 
0.0 for test case 1. The reason for using ao = 0.0 for test case 1 is that there are no 
shock waves present. The sensor is used merely to suppressed high frequency producing 
nonlinear instability associated with the central base scheme in this long time wave 
propagation phenomena. These various numerical methods will be notated as ACM 
or WAV (depending on whether an ACM or wavelet is used as the sensor) with the 
following numbers indicating the order of the spatial interior base scheme for the inviscid 
and viscous terms. For example, ACM66 (WAV66) means the use of sixth-order central 
as the base scheme for both the inviscid and viscous terms, and ACM as sensor (wavelet 
as sensor). These wavelet sensors are computed using the dimension by dimension 
method as discussed in Section 3.5. In order not to introduce additional notation, 
inviscid flow simulations are designated by the same notation, with the viscous terms 


26 



not activated. Computations using the B-spline and the redundant form of Harten’s 
multiresolution method will be notated as “BS” and “RH” as in WAV66-BS and WAV66- 
RH. Computations using entropy splitting are indicated by adding the letters ENT as 

in ACM66-ENT, WAV66-ENT-BS or WAV66-ENT-RH. Computation using S l . : . = 1, 

3 ' 2 

i.e. , the sensor is turned off and the full amount from the upwind TVD dissipation 
portion is used as the filter, will be notated by TVD as in TVD66. The fifth-order 
WENO scheme of Shu [18] will be notated as WEN05. 

For the second and third test cases, we lower the order of the base scheme near 
the boundary points for the boundary scheme. For example, for the sixth-order base 
scheme, we use a fourth-order central scheme two points away from the boundary point 
and second-order central scheme one point away from the boundary point. For the third 
test case, for simplicity, slip wall boundary condition (BC) is used for the lower wall, 
and the upper y- direct ion physical BC is overspecified and nonreflecting BC is not used. 
A uniform Cartesian grid of 80 x 79 is used for test case 1. For test cases 2 and 3, a 
uniform Cartesian grid is used in the .x- direct ion and a mildly stretched Cartesian grid 
is used in the y-direction with the grid size of 101 x 101 and 321 x 81 respectively. 

5.1 Comparison of the Wavelet Sensor with the ACM Sensor 

Before showing the comparison of the high order wavelet based filter scheme with the 
high order ACM based filter scheme, we first show the performance of the wavelet sensor 
using the dimension-by-dimension approach for a 2-D complex flow structure and then 
compare loosely with the ACM sensor. The initial illustration does not involve the 
numerical scheme and only the performance of the wavelet sensor is demonstrated. 
Figure 5.4 shows the computed density and pressure contours at t = 120 by WAV66- 
RH with At = 0.12 for test case 3. Here we only consider these numerical data as a 
given two-dimensional discrete function to be analyzed by the wavelet algorithm. The 
function represents a shock from the upper left corner, impinging on a horizontal shear 
layer in the middle of the domain (See Fig. 5.3). The shock is reflected from the lower 
wall boundary. For more details about the problem, see Yee et al. [24, 25]. 
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Density 



Pressure 



Fig. 5.4. 2-D Testing discrete function, 

(density and pressure contours at t = 120 from test case 3). 


Figure 5.5 shows contours of the estimated Lipschitz exponent a for the function in 
Fig. 5.4. The value a was computed here from three levels (mo = 3) of the wavelet 
algorithm, using the wavelet coefficient 


w m,j,k — 



where the one dimensional coefficients were computed by the multiresolution operators 
(3.24) in each coordinate direction. The coefficients were computed for the pressure. 
The top figure in Fig. 5.5 shows a contours on levels from 0.5 to 0.9. The lower figure 
shows the corresponding sensor, a function which is one for a < 0.5 and zero otherwise. 
The wavelet sensor clearly captures the shock and the shear layer. The low a at the 
upper boundary to the right is probably due to mildly unstable boundary conditions at 
the upper boundary. 
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Alpha contours [0.5 0.9] 



Sensor, contour at 0.5 



Fig. 5.5. Top: a contours 0.5 < a < 0.9; Bottom: sensor contour at a = 0.5. by the RH-wavelet. 

We want to emphasize that Fig. 5.5 shows the sensor when applied to a precomputed 
solution at a fixed time. No dynamic behavior was involved (i.e. , the numerical scheme 
is not part of the analysis). Since the ACM sensor has no one-to-one correspondence 
counterparts of Fig. 5.5, no results are shown for the ACM sensor. Next we show in 
Figs. 5.6 and 5.7 results from actually computing the flow using the respective wavelet 
based (WAV66-RH) and ACM based (ACM66) filter schemes for At = 0.12. Figure 5.6 
shows the wavelet sensor applied to the density and pressure at t = 120 in the x and 
y - direct ions, and the square root of the sum of these quantities in the x and y-directions. 
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Density/Pressure based WAV66-RH x-sensor 


10 

>* 0 
-10 


Density/Pressure based WAV66-RH y-sensor 


10 

>* 0 
-10 


Density/Pressure based WAV66-RH sqrt(x 2 +y 2 ) 


10 

^ 0 
-10 


Fig. 5.6. One contour at a = 0.5 of the sensors used by WAV66-RH 
applied to the density and pressure of test case 3. 

Figure 5.7 shows the corresponding contours using the ACM sensor with k = 0.35. There 
is only one contour level plotted. The level value is in the middle of the range, i.e., at the 
average of the maximum and the minimum sensor values. The wavelet sensor was able to 
extract the full features of the flow structure far better than the ACM sensor. Although 
this is the case, as we can see later, the wavelet sensor exhibits accuracy similar to the 
best tuned ACM sensor. This is due partly to the fact that in actual implementation 
formulas (4.4) - (4.5) are used. We are not making full use of the Lipschitz exponent. 
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Density/Pressure based ACM66 x-sensor 



50 100 150 


x 


Density/Pressure based ACM66 y-sensor 



50 100 150 


x 


Density/Pressure based ACM66 sqrt(x 2 +y 2 ) 



50 100 150 


x 

Fig. 5.7. One contour of the sensor, <Sf +i used by ACM66 
applied to the density and pressure of test case 3. 

Again one cannot compare the ACM sensor and the wavelet sensor directly on these 
figures. The wavelet has more flexibility and choices whereas the ACM sensor only 
compares the strength of gradients (rn = p = 1 in (2.8) and (2.9)) between neighboring 
grid points of a chosen physical quantity (or vector). Perhaps the comparisons would 
be more relevant if we were to compare the ACM sensor using different p and m values 
in (2.8) and (2.9) for the different flow features. This involves additional switching 
parameters and is not pursued here. 
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5.2 Comparison Among TVD66, WEN05, ACM66, and WAV66 

Sample computations using the high order wavelet based filter scheme WAV66 compared 
with the high order ACM based filter scheme ACM6G for test cases 1-3 are shown in 
Figures 5.8-5.11. The accuracy of the two wavelet sensors, B-spline wavelet (WAV66- 
BS) or the redundant form of Harten wavelet (WAV66-RH) for test cases 1-3 (results 
not shown) is very similar and the effect on accuracy of the choice of the physical vector 
(density and/or pressure, characteristic variables, o/. + i, or entropy variables W) to be 

sensed is not pronounced. In all cases, no physical problem-dependent parameter has to 
be tuned. The accuracy compared very well with that of the corresponding best tuned k 
for ACM66 for the individual test cases 1-3. In particular, similar accuracy was sustained 
using the redundant form of Harten wavelet sensor and entropy splitting (WAV66-ENT- 
RH) for long time integrations of the vortex convection problems as ACM66-ENT using 
k = 0.01 and At = 0.01. 

Figure 5.8 shows the comparison among TVD66 (5* , = 1), WEN05, ACM66, 

? '2 

and WAV66-RH for test case 3, illustrating the normalized temperature and pressure 
contours at t = 113.16 with k = 0.35 for the nonlinear fields and k = 0.175 for the 
linear fields for the ACM66. For this set of tuned k’s, the solution obtained is very 
accurate and without visible instability. The solution is comparable with the WAV66- 
RH. For test cases 2 and 3, only 50% of the wavelet sensor is applied to the linear fields 
(i.e. , 50% numerical dissipation). The resolution of the WAV66-RH is more diffusive 
if full strength is applied to the linear fields. Observe that the wavelet sensor was 
able to remove the noise generated on the upper boundary due to the overspecified 
BC. Note that the normalized temperature is the most sensitive value to examine for 
accuracy of the schemes. By examining temperature contours, we note that the vortices 
are more diffusive in the WEN05 computations. There is a minor difference in the two 
simulations. The WEN05 code has a built-in nonreflecting BC on the upper y-direction. 
The WEN05 also requires more operations count than the ACM66 or WAV66. 

The long time wave propagation of the inviscid vortex convection problem simi- 
lar to test case 1 poses a different challenge to the numerical method. For long time 
wave propagation of this nature, non-dissipative or low dissipative schemes usually ex- 
hibit high frequency oscillation nonlinear instability at quite early stages of the wave 
propagation process. Although numerical dissipation can suppress the high frequency 
oscillation, if applied incorrectly, the vortex becomes very diffuse at a longer time in- 
tegration and eventually nonlinear instability sets in. Figures 5.9 and 5.10 show the 
long time wave propagation comparison between ACM66-ENT and WAV66-ENT-RH 
for test case 1. Figure 5.9 shows the density profiles at the centerlines y = 0 and at 
x = 5, cutting through the center of the initial vortex, at 20 spatial period increments. 
The time required for one spatial period is T = 10. 
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Fig. 5.9c. Density at the line y = 0 at T = 20 increment for test case 1, WAV66-ENT-RH. 



Fig. 5.9d. Density at the line x = 5 at T = 20 increment for test case 1, WAV66-ENT-RH. 

The time step and grid spacing are At = 0.01 and 80 x 79. Depending on the time 
and spatial discretizations, the grid size and time step, the vortex, after long time 
integrations, can drift away from the centerline. The amount of drift depends on the 
scheme, the grid size and the time step. For the various methods that we studied in 
[25], the amount of drift can be very severe. For the present two methods with the 
indicated time step and the uniform grid spacings, there is only a very slight drift of 
the vortex after a very long time integration. If the computed vortex drifts away from 
the centerline but still preserves the vortex shape and strength, the centerline density 
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profiles do not convey the full information and can be misleading. We complement 
the comparison with snap shots of density contours at different times up to 130 spatial 
periods in Fig. 5.10. 

Figure 5.10 shows the snap shots of the density contours at different times. The 
vortex is convected 10 spatial periods between each plot. The result from the WAV66- 
ENT-RH method is at least as good as the ACM66-ENT method. Here the base scheme 
is applied to the entropy split form of the inviscid flux derivative, in order to reduce 
effects from non-linear instabilities. The results using the same condition and parame- 
ters, but with no entropy splitting of the inviscid flux derivatives, although very stable, 
exhibit smearing of the vortex and severe vertical and horizontal drifts for both ACM66 
and WAY66-RH. See [25] for the ACM66 and ACM66-ENT comparison. The use of 
entropy splitting in conjunction with ACM66 (ACM66-ENT) or WAV66-RH (WAV66- 
ENT-RH) has preserved a horizontally convecting vortex with great accuracy after long 
time integration of 130 (T = 1300) periods. The results use a uniform and not very fine 
grid. To the authors’ knowledge, highly accurate numerical simulation of this problem 
previously reported in the literature were only carried out up to 10 periods of integra- 
tion. 

We would like to point out that the vertical and horizontal drifting (or rather shift- 
ing) of the vortex away from the centerline y = 0 and/or x = 5 is quite common for 
all schemes beyond 30 periods. Depending on the scheme, the amount of numerical 
dissipation and the time step, drifting can occur as early as 5 periods. We believe that 
the drifting is due largely to the spatial numerical dissipation of the scheme provided a 
highly accurate low phase error time integrator is used. 

Figure 5.11 shows the comparison among TVD66, ACM66, WAY66-RH and WAV66- 
BS for test case 2 with At = 0.1. Here two different versions of the WAV66-RH method 
are examined. The one denoted WAV66-RH in the figure has 50% reduced TVD dissi- 
pation on the linear fields. The resolution of WAV66-RH and WAV66-RHb is slightly 
more accurate than WAV66-BS. The result using WEN05 (not shown) is less accurate 
than ACM66 but more accurate than TVD66. 
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Fig. 5.11. Normalized temperature contours for test case 2 
using TVD66, ACM66, WAV66-RH, WAV66-RHb, and WAV66-BS. 
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6 Concluding Remarks 

Improved adaptive numerical dissipation controls over the Yee et al ACM sensor have 
been constructed. The new sensors with improved detection properties are derived from 
multiresolution wavelet based analysis and require slightly more operations count than 
the ACM sensor. There are a variety of wavelets to choose from, depending on the flow 
feature. 

We considered two types of non-orthogonal wavelet basis functions for our 2-D com- 
pressible Euler and Navier-Stokes numerical experiments. One is similar to the B-spline 
wavelet (Mallat &; Zhong [13]) used by Gerritsen & Olsson [3] for grid adaptation and 
the other is a modification of the multiresolution method of Harten [5] as a redundant 
multiresolution wavelet. The B-spline wavelet sensor requires slightly more arithmetic 
operations and a wider grid stencil than the redundant form of Harten wavelet sensor. 
The final form of the wavelet sensor S l i involves mainly nested difference operators 

and least squares fits. From the numerical experiments, it appears that the RH- wavelet 
sensor exhibits a slightly more accurate result than the BS-wavelet sensor. The proposed 
wavelet sensors, unlike the ACM sensor, can detect most of the distinct flow features, 
including turbulence, leading to an automatic selection of the appropriate distribution 
of numerical dissipation. Since distinct Lipschitz exponent values represent distinct 
flow structures, these wavelet sensors are free of physical problem-dependent arbitrary 
parameters for the three test cases presented. They are also good grid adaptation indi- 
cators [3] when compared to the ones commonly used in practice. Consequently, a new 
dual purpose adaptive method is readily available leading to dynamic numerical dissi- 
pation controls and improved grid adaptation indicators. This dual purpose adaptive 
method can also serve as a stand alone option for other numerical schemes. 

In the future, we will explore the full capability of the multiresolution wavelet prop- 
erty. This will include improved wavelet basis functions and their scaling functions for 
high speed compressible shock-turbulence interaction and numerical combustion, and 
an improved switching function other than the one proposed in Section 4. In other 
words, a better use of the Lipschitz exponent information will be implemented. An- 
other possibility is to use an approximation adapted to the regularity of the solution 
similar to the (h,p) finite element method. We will also consider integrating the wavelet 
sensor into a numerical scheme, by using the Lipschitz exponent a in the same way as 
the switching quantity is used in the high resolution shock-capturing schemes. 
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